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PREFACE 


The  work  we  report  here  has  two  primary  objectives;  (1)  to  investigate  the  effects  of  near¬ 
source  heterogeneity  on  the  radiation  from  explosion  sources,  and  (2)  to  develop  improved 
methods  for  modeling  regional  wave  propagation  in  laterally  varying  media. 

The  first  section  of  this  report  is  a  preprint  of  an  article  submitted  to  Bulletin  of  the 
Seismological  Society  of  America,  and  deals  with  nuclear  explosions  in  and  near  underground 
cavities.  Using  numerical  modeling  with  a  2-D  multiple  multipole  method,  we  study  the 
source  radiation  from  an  explosion  in  a  cavity  including  the  effects  of  scattering  from  a 
nearby  cavity.  We  also  study  the  radiation  from  an  explosion  in  a  regular  arrangement  of 
cavities  and  in  a  particular  ‘L’-shaped  cavity  geometry  representing  a  large  tunnel  with  a 
side-drift.  Our  modeling  results  are  used  to  interpret  seismograms  recorded  during  the  Non- 
Proliferation  Experiment  (NPE),  which  showed  strong  azimuthal  variation  in  waveforms 
and  amplitudes  that  can  be  explained  only  partially  by  heterogeneities  along  the  path  of 
propagation  or  near  the  receiver.  Near  the  NPE  source  there  was  a  large  complex  of  tunnels 
and  cavities.  Synthetic  seismograms  for  a  simplified  NPE  geometry  correlate  nicely  with  the 
observed  S  wave  energy  generated  at  or  near  the  source.  By  including  near-receiver  tunnels 
into  our  model,  we  also  obtain  the  observed  coupling  between  the  radial  and  transverse 
component. 

The  second  section  of  this  report  describes  a  new  technique  we  have  developed  for  com¬ 
puting  synthetic  seismograms  in  laterally  varying  media.  This  technique  parameterizes  the 
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model  on  a  (possibly  irregular)  triangular  grid  and  uses  a  discrete  approximation  to  the 
integral  form  of  the  wave  equation.  An  irregular  grid  has  the  advantage  that  it  can  be 
aligned  with  irregular  interfaces  present  in  a  model,  in  particular  an  irregular  free  surface. 
The  method  is  thus  useful  for  the  incorporation  of  surface  topography  into  a  model,  which  is 
essential  for  the  modeling  of  regional  wave  propagation  in  many  important  geographic  areas, 
e.g.  the  Middle  East.  Another  attractive  feature  of  irregular-grid  modeling  is  the  possibility 
of  locally  adapting  the  grid  spacing  to  the  wave  speed,  saving  memory  and  computer  time. 
We  first  demonstrate  the  improved  accuracy  of  the  method  for  irregular  interfaces  through 
a  comparison  with  finite  differences.  We  then  present  two  models  that  show  the  effect  of  a 
rough  surface  topography  on  regional  seismograms. 


Calculating  Seismic  Radiation  Patterns  for  Explosions 

in  and  Around  Cavities 


Abstract 

Seismic  sources  can  be  located  in  an  underground  complex  of  tunnels  and  cavities.  Examples  of 
such  sources  are  chemical  or  nuclear  devices  as  well  as  mining-related  events.  Cavities  are  very 
strong  scatterers  which  affect  radiation  patterns  and  conversions  between  different  wavetypes. 

We  use  a  2-D  multiple  multipole  method  to  calculate  source  radiation  patterns  and  scattering 
from  cavities.  We  illustrate  the  effects  of  distance  between  source  cavity  and  a  scatterer,  the  effect 
of  a  regular  arrangement  of  cavities  and  the  location  of  the  source  therein,  and  the  effect  of  different 
source  mechanisms  for  a  particular  ‘L’-shaped  geometry  representing  a  large  tunnel  with  a  source 
in  a  side-drift. 

Seismograms  recorded  during  the  Non-Proliferation  Experiment  (NPE)  show  strong  azimuthal 
variations  in  waveforms  and  amplitudes  which  can  be  explained  only  partially  by  heterogeneities 
along  the  path  of  propagation  or  near  the  receiver.  Near  the  NPE  source,  there  was  a  large  complex 
of  tunnels  and  cavities.  Synthetic  seismograms  for  a  simplified  NPE  geometry  correlate  nicely  with 
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data  observed  around  the  arrival  of  S-waves  generated  at  or  near  the  source.  By  including  near¬ 
receiver  tunnels  into  our  model,  we  also  obtain  the  observed  coupling  between  the  radial  and 
transverse  component. 

Introduction 

Discrimination  between  seismic  sources  such  as  earthquakes,  quarry  or  mine  blasts,  and  nuclear 
explosions  becomes  increasingly  difficult  as  the  magnitude  of  the  events  decreases  (Blandford,  1981; 
Li  et  ah,  1996).  The  ability  to  extract  useful  information  about  seismic  sources  from  seismic 
waveform  data  is  limited  by  the  complexity  of  the  seismograms  and  the  ability  to  model  these 
complexities.  For  example,  spall,  tectonic  release  processes,  and  scattering  contribute  strongly  to 
far-held  seismic  wave  signals.  The  classic  explosion  source  model,  which  yields  an  isotropic  P-wave 
radiation  pattern,  cannot  reproduce  observed  waveforms. 

Signal-generated  “noise”  is  particularly  treacherous  in  this  respect,  as  it  typically  has  the  tran¬ 
sient  character  and  frequency  content  which  can  easily  be  misinterpreted  as  due  to  some  property 
of  the  source.  For  example,  changes  in  source  cavity  shape  and  source  location  within  the  cavity 
influence  the  radiation  pattern  of  the  explosion  by  decoupling  (Latter  et  ah,  1961).  Although  many 
studies  have  been  made,  numerical  (e.g.,  Gibson  et  al.,  1996)  and  analytical  results  (e.g.,  Ben- 
Menahem  and  Mikhailov,  1995)  apply  different  assumptions  and  produce  different  results.  Zhao 
and  Harkrider  (1992)  developed  solutions  for  the  displacement  held  created  by  an  explosion  source 
located  off-center  within  a  solid  sphere  embedded  in  an  inhnite  solid  medium,  a  model  of  a  fully 
tamped  explosion.  Here  the  loss  of  symmetry  causes  shear  waves  to  be  generated,  even  with  the 
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spherical  cavity.  Other  studies  have  examined  the  effect  of  varying  the  cavity  shape.  Rial  and 
Moran  (1986),  Glenn  and  Rial  (1987),  and  Gibson  et  al.  (1996)  considered  ellipsoidal  and  cylindri¬ 
cal  cavity  shapes,  and,  applying  a  combination  of  numerical  and  approximate  analytical  methods, 
obtained  estimates  of  far-field  radiation  patterns.  They  suggested  that  the  P-wave  radiation  would 
be  strongest  perpendicular  to  an  ellipsoidal  cavity  and  that  a  fairly  strong  S-wave  would  be  gen¬ 
erated  as  the  cavity  aspect  ratio  (ratio  of  length  to  diameter)  became  large.  Stevens  et  al.  (1991) 
and  Murphy  et  al  (1995)  also  examined  an  ellipsoidal  cavity.  However,  in  contrast  with  the  earlier 
results,  they  concluded  that  at  the  frequencies  of  interest  for  regional  wave  propagation,  the  P-wave 
radiation  pattern  would  be  fairly  isotropic  and  S-wave  amplitudes  comparatively  small. 

In  the  above  studies,  the  effect  of  other  heterogeneities  has  been  mostly  ignored.  However,  App 
et  al.  (1995)  found  that  a  horizontal  discontinuity  near  an  explosion  source  generated  enough  S- 
waves  to  influence  regional  waveforms,  and  thus  cause  such  an  explosion  to  appear  more  earthquake¬ 
like.  For  explosions  at  the  Nevada  Test  Site,  Stump  et  al.  (1994)  noted  that  at  very  closed-in 
distances,  the  transverse  component  of  motion  contained  a  source  signature  which  relates  to  some 
type  of  scattering  mechanism.  An  analysis  of  seismic  data  from  the  Non-Proliferation  Experiment 
(NPE)  showed  that  elastic  wave  scattering  near  the  receiver  had  a  significant  effect  upon  the 
recorded  waveforms  (Johnson,  1995).  Kamm  and  Bos  (1995)  calculated  the  effects  of  the  layered 
structure  on  the  NPE  to  study  the  scattering  along  the  path  of  propagation.  Their  results  indicate 
that  the  velocity  structure  alone  cannot  explain  the  scattered  waves. 

Chemical  or  nuclear  devices  are  commonly  exploded  within  underground  complexes  with  various 
explosion  chambers  as  well  as  shafts  and  tunnels  used  for  access  and  instrumentation  (Denny  and 
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Stull,  1994).  Moreover,  there  exist  other  seismic  sources  which  occur  in  or  near  underground 
complexes  such  as  mines:  e.g.  rockbursts  (Bennett  et  al.,  1995),  collapses  (Pearson  et  al,  1996),  or 
mining  related  explosions  (Reamer  et  al,  1992).  But  cavities  are  extremely  strong  scatterers  which 
efficiently  convert  different  wave-types  into  each  other.  The  scatterers  act  as  secondary  sources 
with  different  radiation  patterns  and  different  distance  dependencies  than  the  primary  waves.  The 
discrimination  between  different  seismic  sources  is  made  much  more  difficult  by  these  secondary 
signals. 

The  objective  of  this  paper  is  to  study  how  the  shape  of  the  source  cavity  and  the  scattering 
by  heterogeneities  in  the  near-field  of  the  source,  such  as  tunnels  and  other  cavities,  affect  the 
radiation  of  seismic  waves.  Characterizing  seismic  radiation  from  such  complex  source  geometries 
requires  an  accurate  and  efficient  computational  method  that  will  work  for  complex  arrangements 
of  cavities  with  a  source  in  any  part  of  a  cavity  complex.  For  this  purpose,  we  developed  the  2-D 
Multiple  MultiPole  (MMP)  method  (Imhof,  1996).  In  the  MMP  method,  scattered  wavefields  are 
expanded  into  a  set  of  basis  functions  which  are  solutions  to  the  wave  equation  in  a  homogeneous 
medium.  To  find  the  yet  unknown  weighting  coefficient  for  each  basis  function,  we  choose  discrete 
matching  points  along  the  boundaries  of  the  scatterers  where  we  satisfy  the  boundary  conditions. 
Each  matching  point  provides  two  equations  corresponding  to  vanishing  tractions  which  build  a 
linear  matrix  system  for  the  unknown  weighting  coefficients.  Typically,  the  chosen  basis  functions 
are  scalar  Helmholtz  potentials  of  multipole  form  for  P-  and  S-waves.  In  contrast  to  more  tradi¬ 
tional  methods,  we  use  multiple  expansion  centers  and  thus  multiple  multipoles  for  each  scatterer. 
However,  these  basis  functions  are  not  orthogonal,  leading  us  to  overdetermine  the  linear  matrix 


system  and  solve  it  in  the  least-squares  sense.  The  method  has  been  described  in  detail  in  a  pre¬ 
vious  publication  (Imhof,  1996),  and  in  Appendix  A  we  outline  the  method  as  necessary  for  the 
present  work.  For  a  more  complete  discussion  on  usage  and  implementation,  the  reader  is  referred 
to  the  previous  paper.  For  some  models,  we  will  convert  the  MMP  solutions  into  moment  tensor 
expansions.  The  relation  between  the  MMP  expansion  and  the  moment  tensor  is  presented  m 

Appendix  B. 

The  paper  is  structured  as  follows.  First,  we  illustrate  the  effect  of  distance  between  an  axially 
symmetric  source-cavity  and  a  scatterer.  We  also  show  how  a  regular  arrangement  of  cavities 
affects  the  radiation  patterns  as  a  function  of  source  location.  Second,  we  exemplify  the  effect 
of  asymmetric  cavity  shapes.  A  small  source  cavity  is  placed  near  a  large  tunnel.  We  show  the 
radiation  patterns  for  only  the  source  cavity,  the  tunnel  in  the  near-field  of  an  isotropic  explosion 
source,  and  combined  models  with  different  source  mechanisms.  Finally,  we  apply  our  findings  to 
the  Non-Proliferation  Experiment.  We  interpret  phases  recorded  around  the  arrival  time  of  S- waves 
by  near-source  scattering  and  explain  the  coupling  between  radial  and  transverse  component  by 
near-receiver  heterogeneities. 

Effect  of  Cavities  Near  the  Source 

As  a  first  example,  we  calculate  how  a  2-D  cylindrical  scatterer  in  the  near-field  of  the  source  affects 
the  radiated  waves.  This  simple  example  serves  two  purposes:  First,  as  a  check  for  the  numerical 
method  we  use.  Second,  the  example  shows  that  a  scatterer  in  the  near-field  of  the  source  has 
a  pronounced  effect  on  the  source  signature.  A  schematic  of  the  example  is  shown  in  Figure  1. 
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Two  cylindrical  cavities  with  a  radius  of  7.5  m  are  embedded  in  a  homogeneous  fullspace  with 
density  p  =  1800 kg/m3,  P-wave  velocity  vp  =  2000 m/s,  and  S-wave  velocity  vs  =  1110 m/s.  One 
cavity  (source)  is  assumed  to  be  instantaneously  pressurized  by  an  explosion  (Patton,  1991).  The 
center-frequency  of  the  radiated  wave  is  10  Hz. 

Radiation  patterns  R(<f>)  as  a  function  of  azimuth  <j>  are  calculated  for  radial  and  transverse 
polarizations: 

Ri{4>)  —  f  {ui(ro,<P,r))2  dr,  (1) 

Jo 

where  the  subscript  i  denotes  radial  or  transverse  component  of  the  total  displacement  u(x,  t)  at 
location  x  =  (r0,<j>).  The  distance  r0  is  750  m.  We  do  not  need  to  decompose  the  seismograms  into 
P-  and  S-waves.  Although  the  P-waves  (S-waves)  will  bleed  into  the  transverse  (radial)  component, 
the  radial  (transverse)  radiation  pattern  is  dominated  by  the  P-waves  (S-waves). 

Due  to  the  circular  cross-section  in  this  2-D  example,  the  source  emits  only  P-waves  where 
the  amplitudes  are  independent  of  angle.  However,  the  second  cavity  placed  nearby  perturbs  the 
isotropic  P-wave  radiation  pattern  and  acts  as  a  secondary  S-wave  source  (Figure  2).  The  farther 
the  two  cavities  are  separated,  the  smaller  the  effect  of  the  scatterer  becomes.  Although  the  P- 
wave  radiation  pattern  is  clearly  anisotropic  for  small  separations  (20  m),  it  is  the  S-wave  radiation 
pattern  which  is  affected  the  most.  While  the  S-wave  displacement  is  nearly  comparable  to  the 
P-wave  displacement  at  20  m,  it  vanishes  quickly  with  increasing  separation.  But  as  the  overlay 
of  the  results  from  all  cavity  locations  shows,  the  originally  isotropic  P-wave  radiation  patterns  is 
only  slightly  perturbed  by  the  presence  of  additional  cavities  in  the  near-field  of  the  source. 

The  rapid  decay  of  the  S-wave  conversion  is  presented  in  Figure  3.  Shown  is  the  ratio  between 
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the  maximal  S-wave  radiation  and  the  isotropic  P-wave  radiation  as  a  function  of  separation  between 
the  two  cavities.  While  the  ratio  is  as  high  as  0.7  for  a  distance  of  20  m,  it  decays  rapidly  until  it 
reaches  a  nearly  constant  value  of  0.07  for  distances  of  300  m  and  more.  Hence  for  multiple  cavities, 
we  expect  only  cavities  within  a  dominant  shear  wavelength  (110  m)  from  the  source  to  influence 
the  S-wave  conversion  and  thus  the  S-wave  radiation  pattern.  The  effects  from  cavities  at  larger 
distances  will  be  overshadowed  by  cavities  near-source  and  near-receiver  heterogeneities. 


Effect  of  a  Row  of  Cavities 

Although  we  concluded  from  the  previous  numerical  experiment  that  the  effect  of  other  cavities 
vanishes  rapidly  with  increasing  distance,  there  still  are  interesting  effects  to  be  observed.  The 
geometry  of  a  complex  of  cavities,  drifts,  and  shafts  might  not  be  random  but  systematic  and 
orderly.  A  likely  arrangement  of  explosion  chambers  would  be  an  array  of  evenly  spaced  cavities. 
The  configuration  we  look  at  is  an  array  of  cavities  with  a  radius  of  7.5  m  and  a  spacing  of  50  m. 
Assuming  an  array  of  9  cavities,  we  will  pressurize  each  of  them  in  a  sequence  (Figure  4).  For  each 
source,  we  will  calculate  the  overall  radiation  pattern  as  a  function  of  position  in  the  sequence. 
The  model  parameters  are  the  same  as  in  the  previous  one:  2-D  model,  density  p  =  1800 kg/m3, 
P-wave  velocity  vp  —  2000  m/s,  S-wave  velocity  vs  =  1110  m/s,  and  center-frequency  10  Hz. 

Figure  5  presents  the  radiation  patterns  for  the  source  chambers  1,  3,  and  5.  The  first  explosion  is 
the  most  asymmetric  configuration,  while  the  fifth  explosion  will  be  in  a  symmetrical  one.  Although 
the  P-wave  radiation  pattern  is  anisotropic  for  each  case,  the  overlay  of  the  three  cases  shows 
that  the  P-wave  pattern  hardly  depends  on  the  choice  of  the  source  chamber  within  the  array. 
Clearly,  the  S-wave  pattern  exhibits  more  change.  While  the  pattern  itself  changes  only  slightly,  the 
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amplitudes  in  the  forward  direction  exhibit  the  greatest  effect.  In  this  direction,  the  conversion  to  S- 
waves  doubles  between  source  chamber  1  and  5.  For  source  chamber  5,  the  ratio  of  S-wave  radiation 
to  P-wave  radiation  reaches  as  high  as  0.6.  The  efficiency  of  the  P  to  S  conversion  as  a  function 
of  source  chamber  is  shown  in  Figure  6.  Once  the  source  chamber  is  within  the  array,  the  ratio 
of  S-wave  radiation  to  P-wave  radiation  stays  close  to  0.6  for  all  source  chambers.  The  exceptions 
are  the  first  and  the  last  source  chamber,  each  of  which  has  only  one  neighboring  cavity.  However, 
even  for  these  extremal  geometries,  the  conversion  to  S-waves  is  much  stronger  than  expected  for 
a  separation  of  50  m  between  source  chamber  and  an  adjacent  cavity.  Comparison  of  Figures  3 
and  6  shows  that  the  ratio  is  two  to  three  times  higher  than  expected  for  this  separation  distance. 
Thus,  in  the  array  configuration,  cavities  at  larger  distances  still  contribute  to  the  conversion  due 
to  multiple  scattering  between  the  cavities. 

Effects  of  Cavity  Shapes  and  Source  Types 

The  two  examples  in  the  previous  section  were  both  for  axially  symmetric  cavities.  Here,  we  will 
address  the  effect  of  asymmetrical  cavities.  The  2-D  geometry  we  will  use  is  shown  in  Figure  7. 
The  situations  we  study  are  sketched  in  Figure  8.  The  models  consist  of  (a)  a  source  in  a  small, 
elongated  cavity  (26  x  10  m),  (b)  an  explosive  line  source  near  a  large  cavity  (60  x  18  m),  and  (c) 
both  cavities  together  representing  a  large  access  tunnel  with  a  source  in  a  perpendicular  side  drift. 
The  distance  between  the  two  cavities  is  10  m.  The  geometry  is  not  plane-symmetric.  The  small 
source  cavity  is  offset  by  15  m  from  the  center  of  the  larger  one.  At  a  distance  of  750  m,  we  calculate 
the  radiation  patterns.  As  the  source  cavity  is  not  axially  symmetric,  it  will  radiate  both  P-  and 
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S-waves.  The  model  parameters  are:  density  p  =  1800  kg/m3,  P-wave  velocity  vp  —  2000  m/s, 
S-wave  velocity  vs  =  1110  m/s,  and  center-frequency  10  Hz. 

The  first  numerical  experiment  is  simply  to  determine  the  radiation  pattern  for  the  source  cavity 
only.  Thus,  we  omit  the  large  cavity  for  this  calculation.  As  a  source,  we  choose  an  instantaneous 
pressurization  of  the  cavity.  Figure  8(a)  shows  a  sketch  of  the  model  and  the  resulting  radiation 
patterns  of  P-  and  S-waves.  The  ratio  between  peak  amplitudes  of  P-  and  S-waves  is  very  close 
to  unity.  As  expected,  we  obtain  a  large  S-wave  component  with  a  radiation  pattern  typical  for 
a  double-coupled  system  of  forces.  For  the  P-waves,  we  see  a  combination  of  isotropic  radiation 
(point  source)  and  quadrupole  radiation  (double-coupled  force  system).  Both  in  our  calculation 
and  those  of  Rial  and  Moran  (1986),  Glenn  and  Rial  (1987),  or  Gibson  et  al.  (1996),  the  P-wave  ra¬ 
diation  is  strongest  perpendicular  to  the  elongated  cavity-axis.  In  contrast  to  these  results,  Stevens 
et  al.  (1991)  and  Murphy  et  al.  (1995)  stated  that  at  the  frequencies  of  interest  for  regional  wave 
propagation,  the  P-wave  radiation  patterns  would  be  fairly  isotropic  and  the  S-wave  amplitudes 
would  be  comparably  small  based  on  their  finite-stress  amplitude  calculations.  For  small  sources 
placed  in  cavities,  elastic  behavior  would  be  a  valid  approximation,  and  our  results  should  represent 
the  radiation  patterns. 

As  a  next  step,  we  determine  the  effects  of  the  large  access  tunnel.  We  replace  the  source  cavity 
with  an  isotropic,  explosive  line  source  located  at  the  center  of  gravity  of  the  source  cavity.  With 
this  arrangement,  there  is  no  interaction  between  the  scatterer  and  the  source.  A  schematic  of 
the  model  and  the  resulting  radiation  patterns  are  shown  in  Figure  8(b).  The  P-wave  radiation 
is  strongly  affected  by  the  tunnel.  The  tunnel  effectively  blocks  the  propagation  of  P-waves.  The 
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radiation  pattern  exhibits  some  degree  of  symmetry.  However,  the  patterns  are  not  aligned  with 
the  scatterer  due  to  the  asymmetrical  configuration  of  source  and  scatterer. 

The  S-wave  radiation  pattern  still  resembles  a  deformed  double-couple.  Interestingly,  the  ratio 
of  S-  to  P-waves  remains  very  close  to  unity.  These  two  examples  clearly  show  that  both  the 
geometry  of  the  source  cavity  as  well  as  additional  cavities  in  the  near-field  of  the  source  generate 
similar  effects,  i.e.  conversion  of  P-  to  S-waves  and  angularly  asymmetric  radiation  patterns. 

We  now  calculate  the  combined  model  where  we  have  an  instantaneously  pressurized  source 
chamber  with  a  large  scatterer  nearby.  Thus,  the  source  chamber  radiates  P-  and  S-waves  as 
shown  in  Figure  8(a).  However,  these  fields  interact  now  with  the  large  access  tunnel.  We  expect 
the  symmetry  in  the  radiation  patterns  to  diminish  which  is  exactly  what  happens.  Figure  8(c) 
shows  the  resulting  radiation  patterns  for  the  combined  model.  The  P-wave  radiation  pattern  still 
resembles  the  original  pattern  as  radiated  from  the  source,  although  there  is  less  P-wave  radiation. 
Also,  the  symmetry  of  the  radiation  pattern  is  broken.  The  big  difference  lies  in  the  S-wave 
radiation.  The  S-wave  radiation  is  enhanced  by  a  factor  of  at  least  1.5.  There  is  a  significant 
increase  in  S-wave  radiation  in  the  diagonal  direction.  To  examine  the  radiation  pattern  further, 
we  convert  the  MMP  solution  into  a  moment  tensor.  The  relation  between  MMP  solution  and 
moment  tensor  is  derived  in  Appendix  B.  This  conversion  will  allow  us  to  relate  our  numerical 
results  to  observations  where  the  moment  tensor  is  routinely  estimated  (e.g.  Patton,  1991;  Johnson 
and  McEvilly,  1994). 

The  moment  tensor  expansion  describes  a  complex  seismic  source  as  a  linear  superposition  of 
an  explosion  source  p,  two  point  forces  and  two  double-couple  force-systems  m,  n.  All 
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components  are  located  at  the  origin.  Customarily,  higher  order  source  components  are  neglected. 
The  components  can  either  be  represented  as  a  function  of  frequency  or  as  a  function  of  time.  In 
the  latter  case,  we  obtain  the  source-time-history.  However,  there  also  exists  a  complementary  set 
of  components.  The  first  component  is  the  rotational  source  s  corresponding  to  a  torque.  The 
remaining  ones  are  non-physical  forces  /i,  /2,  and  non-physical  double-couples  rh,  n.  Their  P-  and 
S-wave  radiation  patterns  are  aligned  instead  of  interlaced  as  for  physical  ones.  These  patterns 
violate  assumptions  underlying  linear,  elastic  wave  theory.  We  choose  the  Haskell’s  source-time 
function  with  B  =  —0.1  and  k  =  25.0s_1  (Haskell,  1967). 

The  different  components  of  the  moment  tensor  expansion  are  shown  in  Figure  10(a).  The 
moment  tensor  is  strongly  deviatoric  which  corresponds  to  a  large  contribution  of  the  m  double¬ 
couple.  Interestingly,  the  presence  of  the  large  scatterer  in  the  near-field  of  the  source  induces  a 
torque  as  indicated  by  the  s  component.  Numerical  experiments  showed  that  moving  the  scatterer 
downwards  into  a  more  symmetric  position  diminishes  the  torque  component  s.  Similarly,  moving 
the  scatterer  into  the  far-field  or  reducing  its  size  decreases  its  strength.  Thus,  the  presence  of  the 
torque  term  is  caused  by  the  asymmetric  geometry.  Both  the  isotropic  component  p  as  well  as  the 
off-diagonal  double-couple  n  contribute  very  little  compared  to  the  deviatoric  double-couple  m.  As 
expected,  the  physical  point  forces  /i,  /2  and  the  non-physical  components  /i,  /2,  hi ,  and  h  do  not 
contribute. 

An  interesting  extension  of  the  previous  discussion  is  to  change  the  source  mechanism.  So  far, 
we  only  used  isotropic  explosion  sources  or  instantaneous  pressurization  of  cavities.  Neither  is  a 
realistic  model  for  underground  mining  where  explosive  charges  commonly  are  placed  in  drill-holes 
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at  the  front  of  drifts.  To  study  the  difference,  we  prescribe  tractions  only  along  one  front  of  the 
source  chamber.  This  source  type  corresponds  to  a  point  force  perpendicular  to  the  front.  Figure  9 
presents  a  sketch  of  the  model  and  the  resulting  radiation  patterns.  Clearly,  the  point  force  is  the 
main  contribution  for  both  P-  and  S-wave  radiation.  However,  the  patterns  are  asymmetric  due  to 
interaction  of  the  waves  with  the  large  scatterer.  These  findings  are  consistent  with  the  moment 
tensor  shown  in  Figure  10(b)  where  the  largest  contribution  comes  from  the  point  force  f\.  Again, 
the  non-physical  terms  fi,  f 2,  and  n  vanish. 

The  radiation  patterns  for  the  two  different  source  mechanisms  ‘instantaneous  pressurization’ 
and  ‘frontal  explosion’  are  clearly  very  different.  The  presence  of  the  access  tunnel  (large  scatterer) 
makes  the  radiation  patterns  deviate  dramatically  from  the  ideal  ones. 


Application  to  the  Non-Proliferation  Experiment 

On  September  22,  1993,  the  Non-Proliferation  Experiment  (Denny  and  Stull,  1994)  was  conducted 
at  the  Nevada  Test  Site  (NTS).  A  large  chemical  charge  equivalent  to  1.1  kt  TNT  was  detonated 
in  a  cavity  off  from  the  N-Tunnel  in  Area  12.  The  purpose  of  the  experiment  was  to  quantify  the 
effects  of  a  chemical  explosion  of  magnitude  comparable  to  that  of  nuclear  explosions.  The  charge 
was  emplaced  in  a  cylindrical  cavity  of  nominal  height  of  5.2  m  and  diameter  15.2  m  at  a  nominal 
depth  of  389  m.  The  source  was  embedded  in  zeolitized  tuff  with  vp  =  2400.0  m/s,  vs  =  1180.0  m/s, 
and  p  =  1800.0 kg/m3  (Kamm  and  Bos,  1995).  The  local  stratigraphy  was  determined  from  three 
vertical  drill  holes  in  a  range  of  300  m  from  the  source. 

Figure  11  shows  seven  seismograms  recorded  at  the  surface  within  a  range  of  558  —  920  m  from 
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the  center  of  the  NPE  explosion  at  different  azimuths.  Range  and  azimuth  are  defined  in  Figure  11 
and  Table  I  (Kamm  and  Bos,  1995).  Clearly,  the  first  motions,  which  correspond  to  the  primary 
P-wave,  are  along  the  radial,  upward  direction.  The  remainder  of  the  seismograms  is  composed  of 
reflections  and  conversions  of  the  initial  wave.  Olsen  and  Peratt  (1994)  speculated  that  scattering 
processes  due  to  medium  heterogeneities  cause  this  effect.  Johnson  (1995)  explained  these  events 
with  scattering  from  heterogeneities  near  the  receivers.  App  and  Brunish  (1991)  explained  similar 
variations  by  reflections  from  a  dipping  low  velocity  layer.  Although  there  is  a  slight  dip  present  in 
the  local  strata  (Baldwin  et  a/.,  1994),  we  do  not  expect  azimuthal  variations  of  such  a  strength. 
We  believe  that  another  possible  scattering  mechanism  -  other  cavities  close  to  the  source  -  may 
also  have  contributed.  For  obvious  reasons,  the  source  cannot  be  placed  in  pristine  material.  The 
source  cavity  will  necessarily  be  surrounded  by  a  large,  complex  network  of  caverns,  drifts,  and 
shafts.  These  cavities  can  be  of  dimensions  comparable  to  the  dominant  wavelength.  Furthermore, 
some  cavities  are  located  within  a  dominant  wavelength  from  the  source.  Hence,  it  is  important 
to  recognize  these  caverns  as  strong  scatterers,  acting  as  secondary  sources  contributing  to  the 
received  signals. 

We  used  the  2-D  MMP  technique  (Imhof,  1996)  to  analyze  the  effect  of  the  caverns  and  tunnels. 
The  complicated  NPE  geometry  shown  in  Denny  and  Stull  (1994)  was  simplified  to  the  model  shown 
in  Figure  12.  Due  to  the  2-D  geometry,  each  cavity  extends  to  infinity  along  the  z-direction.  To 
calculate  seismograms,  we  choose  a  source- time  function  as  suggested  by  Haskell  (1967)  for  an 
underground  explosion  where  B  =  -0.075  and  k  =  30.0  s-1  which  yields  a  source-time  function 
comparable  to  that  observed  by  Johnson  and  McEvilly  (1994).  Figure  13  shows  the  synthetic 
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moment  tensor  as  a  function  of  time.  The  moment  tensor  is  obtained  from  the  MMP  solution  by 
applying  the  conversion  given  in  Appendix  B.  The  synthetic  moment-tensor  and  the  observed  one 
compare  well.  Both  contain  only  small  non-isotropic  components. 

At  a  distance  of  750  m,  we  calculate  the  radial  and  transverse  accelerations  as  a  function  of  time 
and  azimuth.  Figure  14  shows  the  synthetic  seismograms.  There  is  very  little  bleeding  of  energy 
from  the  radial  component  to  the  transverse  one  and  vice  versa.  This  suggests  that  the  radial 
component  contains  mainly  the  P-wave,  while  the  transverse  one  contains  the  S-wave.  Clearly,  for 
waves  propagating  toward  the  east  (Figure  12),  the  pulse  broadens  due  to  the  reflection  off  the 
cavities.  In  contrast,  P-waves  propagating  westward  arrive  with  a  delay  caused  by  the  obstruction 
from  the  tunnels.  For  S-waves,  we  obtain  a  different  picture.  S-waves  propagating  east  arrive  later 
due  to  their  increased  path  from  the  scattering  tunnels.  Waves  propagating  toward  the  northwest 
arrive  earlier  due  to  their  shorter  path.  Unfortunately,  the  dependence  of  amplitude  on  azimuth  is 
not  easily  seen  on  the  seismograms. 

Figure  15  shows  the  corresponding  radiation  patterns  as  defined  in  equation  (1).  The  major 
part  of  the  P-wave  energy  radiates  toward  the  unobstructed  southeast.  S-waves  appear  to  prop¬ 
agate  predominantly  in  the  northeast  and  southward  directions.  Also,  there  are  two  lobes  which 
correspond  to  S-waves  generated  at  the  south  end  of  the  access  drift  and  in  the  gaps  between  the 
three  cavities. 

Figure  16  shows  a  comparison  between  observed  accelerations  recorded  in  the  tunnel  complex 
and  the  synthetics  neglecting  the  effect  of  near-receiver  heterogeneities.  The  observed  accelerations 
are  plotted  in  [m/s2].  For  the  comparison,  the  observed  traces  are  bandpass  filtered  to  contain  a 
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frequency  content  similar  to  the  synthetics.  As  a  filter,  we  use  the  Fourier  transform  of  a  zero-phase 
Ricker  wavelet  with  the  spectrum  5(/)  =  (///o)2  exp{-(///o)2}  where  /o  =  27Hz.  All  synthetic 
traces  are  scaled  equally,  preserving  the  true  amplitude  ratios  between  the  different  receivers  and 
components. 

For  all  four  stations,  synthetics  and  observations  correlate  well  for  early  arrivals  on  the  radial 
component.  These  arrivals  correspond  to  the  initial  P-wave  and  its  reflections  from  heterogeneities 
along  the  path  of  propagation.  On  the  transverse  component,  observed  and  synthetic  traces  corre¬ 
late  well  around  the  arrival  of  S-waves  generated  near  the  source.  These  events  propagated  most 
of  their  path  with  S-wave  velocities.  Not  only  are  the  envelopes  of  the  amplitudes  comparable,  but 
also  their  phases.  In  the  synthetics,  these  S-waves  dominate  the  transverse  component. 

However,  these  synthetics  do  not  explain  the  prominent  transverse  arrivals  at  or  near  the  P- 
wave  arrival  time.  The  accelerometers  were  grouted  into  small  holes  near  access  tunnels  which  act 
as  near-receiver  scatterers.  To  show  their  effects,  we  repeat  the  calculations  but  including  near¬ 
receiver  tunnels.  As  shown  in  Figure  17,  synthetics  and  observations  correlate  well  for  early  arrivals 
on  both  components  for  all  stations.  The  applied  scaling  is  the  same  as  in  Figure  16.  Although 
these  arrivals  propagated  predominantly  as  P-waves  on  the  radial  component,  the  near-receiver 
heterogeneity  couples  these  events  to  the  transverse  component  and  thus  to  S-waves.  Similarly, 
transverse  arrivals  are  imprinted  onto  the  radial  component.  In  conclusion,  both  components 
contain  large  S-wave  contributions  caused  by  near-source  scattering.  The  components  are  coupled 
by  near-receiver  scattering. 
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Summary 


We  have  investigated  the  effects  of  scattering  by  cavities  near  the  source  or  near  the  receivers.  Our 
numerical  analysis  showed  that  scattering  in  the  near-field  of  the  source  has  a  significant  effect 
upon  both  the  observed  waveforms  and  the  radiation  patterns.  In  our  examples,  the  ratio  between 
peak  amplitudes  of  P-  and  S-waves  ranges  from  10  to  70%  depending  on  distance  and  arrangement 
of  source  and  additional  cavities. 

Seismograms  recorded  during  the  Non-Proliferation  Experiment  exhibit  a  strong  variation  in 
their  signals  as  a  function  of  receiver  location.  We  calculated  some  of  these  traces  using  a  simplified 
2-D  model  of  the  underground  test  site.  We  found  that  the  transverse  component  matches  the 
observed  data  fairly  well  around  the  arrival  time  of  S-waves  generated  by  cavities  in  the  near-field 
of  the  source.  However,  early  arrivals  on  the  transverse  component  are  due  to  the  scattering  of 
P-waves  by  near-receiver  heterogeneities.  Similarly,  S-waves  generated  near  the  source  scatter  from 
the  transverse  component  into  the  radial  one. 

We  believe  that  near-source  heterogeneities  have  a  strong  effect  on  the  radiation  patterns  and 
waveforms  received  at  larger  distances.  However,  these  results  were  obtained  from  a  2-D  model.  In 
2-D,  an  extended  cavity  acts  as  a  barrier.  In  3-D,  the  cavity  is  just  a  scatterer.  Hence,  we  expect 
some  differences  between  the  2-D  model  and  the  3-D  experiment.  In  the  future,  we  plan  to  expand 
the  MMP  method  from  2-D  to  3-D  which  will  allow  us  to  address  these  issues  more  rigorously. 
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Appendix  A  Elastic  Multiple  Multipole  Expansions 

The  calculations  in  this  paper  are  done  using  a  multiple  multipole  expansion  method  (Imhof,  1996). 
This  appendix  presents  an  introduction  to  the  method.  The  medium  is  assumed  to  be  a  homoge¬ 
neous  fullspace  with  D  embedded  cavities.  In  the  frequency  domain,  the  displacement  u(x,w)e 
of  an  elastic  P-SV  wave  travelling  in  the  two-dimensional,  homogeneous  region  is  described  by  (Pao 

and  Mow,  1973) 

-i-VV  •  u  —  iv  xVxu  +  u  =  0  (A-l) 

k 2 

where  we  suppressed  the  harmonic  time  factor  eiwt.  We  also  introduced  the  wave  numbers  k  =  u/vp 
and  l  -oj/vs  for  a  particular  frequency  u,  P-wave  velocity  vp,  and  S-wave  velocity  vs. 

Instead  of  using  the  displacement  u(x,  u)  directly,  we  separate  it  into  parts  (Pao  and  Mow, 

1973;  Imhof,  1996): 


u(x,u;)  =  u$(x,cj)  +  u^(x,w)  + uinc(x,o;)  +  ue(x,u;).  (A-2) 

The  first  term  denotes  the  P-waves  generated  by  the  cavities.  Similarly,  the  second  term  denotes 
the  S-waves  radiated  by  the  cavities.  The  third  term  urac(x,u;)  is  used  for  sources  outside  the 
cavities,  or  incident  wavefields  which  are  known.  Because  and  u*  depend  on  the  geometry  of 
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the  cavities  and  the  source  mechanism,  we  will  expand  them  into  multiple  multipole  series.  We  will 
need  to  truncate  the  expansions  after  a  finite  number  of  terms.  Therefore,  we  also  add  an  error 
term  u£  to  (A-2). 

The  displacement  fields  u$  and  u*  are  expanded  as  follows: 

D  Pd  +N 

u*(x,w)  =  ^  j  ^  dpnd  u^d(x,  xpd,  k,  u) ,  (A-3a) 

d=  1  p— 1  n=—N 
D  Pd  +N 

u'i'(x,a;)  =  ^2  bpndUpnd(x,xpd,l,oj) ,  (A-3b) 

d= 1  p= 1  n~—N 

where  apnd  and  bpnd  are  yet  unknown  weighting  coefficients.  The  expansion  functions  u*nd  and 
u^d  are  solutions  to  (A-l).  We  have  an  expansion  for  each  cavity  1  <  d  <  D.  In  each  cavity  d, 
we  place  Pd  different  multipoles  of  orders  —  N  to  +N  centered  at  Each  summation  over  n 
corresponds  to  a  multipole.  Because  we  have  more  than  one  multipole  per  cavity,  the  scheme  is 
named  multiple  multipoles  expansion.  The  expansion  functions  are  defined  as: 

ujL i(x,  Xpd ,  k,  u)  =  ( k  |x  -  xpd|)  em^d  (A-4a) 

u pnd(x,  xpd ,  Z,  o>)  =  V  x  ytfjj  (Z  |x  -  Xpd|)  (A-4b) 

Each  expansion  function  of  (A-4)  satisfies  the  wave  equation  (A-l).  For  each  function,  e.g.  u^nd, 
we  obtain  the  corresponding  stress  tensor  T pnd: 

T %d  =  \C  :  (VuJL  +  u*ndV) ,  (A-5) 

where  C  is  a  fourth-order  tensor  containing  the  elastic  parameters  (Ben-Menahem  and  Singh,  1981). 
On  the  boundary  of  the  cavities,  the  total  traction  t(x,w)  =  T(x,w)  •  n(x)  has  to  vanish.  Thus, 
we  find  a  linear  set  of  equations  by  satisfying  the  boundary  condition  on  discrete  matching  points 
x  at  the  boundary  of  the  cavities.  Along  each  cavity  d ,  we  have: 
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D  Pd  +N  f  'l 

t‘“(x,W)  +  E£  E  <W‘ 5rf(*,w)  +  «Mt^(x,u.)|+f  =  tr(x,«)  (A-6) 

d=l  p=l  n—~~N  ' 

where  we  included  a  prescribed  traction  tfc(x,u;)  along  the  boundary  of  cavity  d.  We  can  use  this 
to  include  instantaneous  pressurization  sources  into  our  models.  Thus,  sources  can  be  represented 
either  by  an  incident  wavefield  uinc(x,u;)  and  its  conjugate  tractions  tinc(x,u/)  or  by  prescribing  the 
traction  t ^(x,^)  along  the  cavity  d.  For  simple  models,  either  uznc(x,  w)  or  t^rc(x,oj)  vanishes. 
Generally,  the  resulting  system  of  linear  equation  needs  to  be  solved  in  the  least-squares  sense  by 
choosing  more  matching  points  than  expansion  functions. 

Appendix  B  Moment  Tensor  from  Multiple  Multipole  Expansions 

MMP-expansions  of  the  form  (A-3)  are  useful  for  numerical  calculations.  However,  we  can  also 
translate  them  into  a  more  traditional  moment  tensor  expansion  which  is  routinely  calculated  from 
observed  data  (e.g.  Patton,  1991).  Graf’s  theorem  (also  called  addition  theorem:  Abramowitz  and 
Stegun,  1964,  Eqn.  9.1.79)  allows  us  to  move  an  expansion  function  to  the  origin.  The  different 
angles  and  vectors  are  defined  in  Figure  18. 

+oo 

H{n\(k  jx  -  xpd|)  ein^d  =  s(n)  s(n  -  m)  s(m)  J\n-m\(k  |xpd|)  e_l(n_m)9 *  Hjm|(k  lxl)  e*m 

TTL—  —  OO 

(B-l) 

Because  we  chose  a  different  normalization  of  the  Hankel  functions,  we  need  to  include  additional 
sign  functions  s(n)  defined  as  s{n)  =  1  for  n  >  0  and  s(n)  =  (-l)n  for  n  <  0.  J\n\(kr)  denotes 
the  Bessel  function  of  order  n.  The  above  version  of  the  translation  theorem  is  only  valid  for 
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|x|  >  |xpd|.  Thus,  we  move  the  complete  right-hand  side  of  (A-3)  to  the  origin  which  collapses  the 
MMP  expansions  into  simple  multipole  (SMP)  expansions  valid  in  the  far-field. 

D  Pd  +N  (  } 

E  E  E  |  apndVH\n\(k  |x  -  Xpd|)  ein^pd  +  bpndV  x  yHln](l  |x  -  xpd\)  e"***  j  = 

d=l  P—1  TL  —  —  N  ^ 

2  {omVfrH(fc|x|)e^  +  /3UVxyfrw(Mx|)c^|  (B-2) 

m— — oc  ' 

In  the  far-field,  the  wavefield  can  also  be  expressed  by  a  moment  tensor  (MT)  expansion  (Backus 
and  Mulcahy,  1976a, b)  describing  the  source.  However,  the  MT  expansion  corresponds  only  to 
one-half  of  the  SMP  expansion  (B-2).  The  other  one-half  of  (B-2)  correspond  to  torques  and 
non-physical  force-systems.  They  are  non-physical  because  they  violate  the  underlying  linear  wave 
theory.  Their  P-  and  S-radiation  patterns  are  aligned.  In  contrast,  the  patterns  are  interlaced 
for  physical  force  systems.  The  Green’s  dyadic  G(x,u>)  for  a  source  at  the  origin  is  defined  to  be 
(Ben-Menahem  and  Singh,  1981) 

G(x,  w)  =  ^  j/2It/>(x)  +  VV  fy(x)  -  <Kx)]  }  ,  (B-3) 

where  I  is  the  identity  tensor.  <£(x)  and  ip(x)  denote  the  scalar  Green’s  functions: 


<^(x)  =  —Ho(k  |x|), 


iP(x)  =  -H0(l  |x|). 


(B-4) 


By  orthogonality,  we  define  another  dyadic  G(x,  u>): 


G(x,  oj)  =  — ~  I  •  VV^(x)  }  where  I  = 
pur  1 


>} 


^0 


V1  V 


(B-5) 
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Using  G(x,u>)  and  G(x,w),  we  write  the  wavefield  u(x,u>)  as: 

u(x,o>)  =  VG(x,u)  •  Ip(w)  +  VG(x,«)  •  is(w)  + 

G(x,w)  •  f(w)  +  G(x,u>)  •  f(w)  + 

VG(x,  u)  •  M(u>)  4-  VG(x,  lj)  ■  M (u)  +  . . .  (B-6) 


p  corresponds  to  an  explosion  source,  f  is  a  point  force.  Each  quantity  with  a  tilde  is  either  a  non¬ 
observed  source-type  (s  corresponds  to  a  torque),  or  are  non-physical,  mathematical  constructs 
without  application  in  seismology.  Both  moment  tensors  M(o>)  and  M(u;)  are  symmetric  and 


deviatoric: 


M(w)  = 


m(u)  n[u) 
n(uj)  -m(u) 


M(u)  = 


m(u ))  h(uj) 

h(u) 


(B-7) 


m  and  n  correspond  to  two  double-couple  systems  which  are  rotated  by  45  0  with  respect  to  each 
other.  The  evaluation  of  (B-6)  yields  that  the  first  two  terms  depend  only  on  Hankel  functions  of 
zeroth  order.  The  next  two  terms  depend  on  Hankel  functions  of  the  first  order,  and  the  final  two 
terms  depend  on  Hankel  functions  of  the  second  order.  Thus,  we  can  equate  (B-2)  and  (B-6)  which 
allows  us  to  find  linear  relations  between  (oo, Pq)  and  (p,s)>  between  (<*±i, P±i)  and  (/ij/2>/i>/2)) 
and  between  (a±2,0±2)  and  (m,  n,  m,  n). 


-4 


(B-8a) 

(B-8b) 
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fi  =  (2  +  2i)  {k(iP+ 1  +  P-i)  ~  K^a+i  +  a-i)} 

(B-9a) 

2 

h  —  (2  +  2i)^-{/(a+i  +  -  k(/?+i  + 

(B-9b) 

2 

/i  =  (2  +  2i)^-{^(a+i  +  iot-i )  +  k(i/3+ 1  +  P-i)} 

(B-9c) 

2 

h  =  (2  +  +  a-i)  _  HP+ 1  + 

(B-9d) 

2 

,  =  -(2  +  2i)pj2  {fc2(^+2  +  ^-2)  -  *2(*<*+2  +  0-2)} 

(B-lOa) 

,  =  -(2  +  2?)  {{2(a+2  +  *a-2)  ~  *  (P+2  +  ip- 2)} 

(B-lOb) 

,  =  -(2  +  2i)  {l2(a+2  +  iot-2 )  +  k  (ip+2  +  P-2) } 

(B-lOc) 

,  —  —(2  +  2x)^2]2  {^2(^a+2  +  <^-2)  “  k2(fi+2  +  iP-2)} 

(B-lOd) 

Once  a  solution  to  a  problem  is  obtained  by  solving  the  linear  system  for  the  weighting  coefficients 
apnd  and  bpnd,  it  can  be  converted  to  a  moment  tensor  using  relations  (B-8),  (B-9),  and  (B-10).  The 
moment  tensor  can  directly  be  used  to  show  the  frequency  dependence  of  the  different  components. 
Alternatively,  the  results  can  be  Fourier  transformed  to  yield  the  source-time  history  as  presented 
in  Figures  10  and  13. 
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Table  I:  NPE  accelerometer  stations  used  in  this  paper.  Listed  are  the  gauge  ID,  the  radial  range 
from  the  source,  the  source-to-station  azimuth  (degrees  East  of  North),  the  fielding  organization, 
and  an  indication  whether  the  gauge  was  located  at  the  surface  or  in  a  cavity  at  the  same  depth 
as  source  measuring  the  freefield  (Kamm  and  Bos,  1995). 
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Figure  1:  Geometry  for  the  first  numerical  experiment:  a  cylindrical  cavity  is  instantaneously 
pressurized  by  an  explosion.  At  a  certain  distance  in  the  near-field,  there  is  a  second  cavity 
scattering  and  converting  wave-modes.  The  total  radiation  pattern  is  then  measured  at  a  distance 
of  750  m  from  the  source. 
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Figure  4:  Array  of  9  cavities  with  a  spacing  of  50  m.  Each  cavity  is  used  once  as  a  source  chamber. 
The  case  shown  is  the  third  one  in  the  sequence.  The  total  radiation  pattern  is  measured  at  a 
distance  of  750  m  from  the  source  chamber. 
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Figure  9:  Radiation  patterns  for  P-waves  (dashed)  and  S-waves  (solid):  The  geometry  is  the  same 
as  in  Figure  8(c).  But  the  seismic  source  is  an  explosion  at  the  front  of  the  source  cavity  as  in  an 
underground  mining  or  construction  environment. 
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Figure  10:  Components  of  the  moment  tensor  expansion  as  a  function  of  time  for  (a)  the  instan¬ 
taneous  pressurization  model  shown  in  Figure  8(c),  and  (b)  the  frontal  explosion  model  defined  in 
Figure  9.  The  components  of  the  moment  tensor  are  explained  in  Appendix  B. 


Figure  11:  Accelerations  recorded  at  a  range  of  558  —  920  m  from  the  center  of  the  NPE  explosion. 
Shown  are  the  vertical  (solid),  radial  (dashed),  and  transverse  component  (dotted).  The  data  was 
acquired  by  LLNL  and  LANL  (Denny  and  Stull,  1994). 
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Figure  14:  The  radial  and  transverse  component  of  the  acceleration  as  calculated  for  the  model 
shown  in  Figure  12. 


Simulation  of  Elastic  Wave  Propagation  in  Models 
Containing  Irregular  Interfaces  Parameterized  in 

Irregular  Grids 


Abstract 

This  study  demonstrates  the  advantages  of  a  recently-developed  irregular-grid  modeling  technique 
(Nolte  et  al.,  1996).  This  technique  can  model  irregular  interfaces  more  accurately  than  a  standard 
regular-grid  finite-difference  method.  We  show  this  by  comparison  of  both  methods  for  a  simple 
model  containing  a  sloping  interface.  While  the  discrete  approximation  to  the  sloping  interface 
results  in  numerical  inaccuracies  for  the  finite-difference  method,  the  irregular-grid  technique  pro¬ 
duces  superior  results.  We  then  show  that  the  method  can  also  be  applied  to  a  free  surface  with 
irregular  topography,  suggesting  that  it  may  be  a  valuable  alternative  to  existing  finite-difference 
free-surface  algorithms.  Finally,  we  demonstrate  the  effect  of  irregular  surface  topography  on  re¬ 
gional  wave  propagation  with  a  Middle  East  example. 
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Introduction 


It  is  well  known  that  finite-difference  techniques  (e.g.,  Alford  et  al.,  1974;  Virieux,  1984,  1986)  suffer 
from  inaccuracies  in  modeling  irregular  interfaces,  due  to  their  discrete  approximation  to  the  model. 
For  example,  a  dipping  interface  has  to  be  approximated  by  discrete  steps.  This  approximation 
can  result  in  diffractions  from  these  “stair  steps”  (Muir  et  al.,  1992).  While  these  diffractions  can 
often  be  eliminated  by  an  isotropic  interpolation  of  the  medium  parameters  across  the  interface, 
more  severe  problems  occur  when  larger  impedance  contrasts  are  present  at  an  interface.  Muir 
et  al.  (1992)  presented  an  effective-medium  approach  to  model  irregular  interfaces.  Even  though 
there  is  no  rigorous  theoretical  justification  for  their  approach,  their  numerical  results  show  that  it 
works  well  in  cases  where  isotropic  interpolation  fails.  Their  method  is  limited  to  interfaces  with 
solid  material  on  either  side,  and  thus  cannot  be  applied  to  an  irregular  free  surface. 

Incorporating  surface  topography  into  finite-difference  modeling  codes  has  been  the  subject  of 
various  previous  investigations  (e.g.,  Jih  et  al.  1988;  Hestholm  and  Ruud,  1994;  Robertsson,  1996). 
As  shown  by  Robertsson  (1995),  various  methods  have  different  shortcomings.  For  example,  grid- 
distortion  techniques  (Tessmer  et  al.,  1992;  Hestholm  and  Ruud,  1994;  Hestholm,  1996)  become 
unstable  at  steep  slopes.  Imaging  techniques  (Levander,  1988;  Robertsson,  1996)  appear  to  be 
more  flexible  in  handling  irregular  topography.  However,  Robertsson  (1996)  was  forced  to  use  a 
very  fine  gridding  in  order  to  avoid  the  stair-step  problems  at  the  free  surface.  Thus,  an  accurate 
representation  of  an  irregular  free  surface  is  computationally  expensive,  and  alternative  modeling 
schemes  are  desirable.  One  approach  could  be  rectangular  variable  grid  methods  (Jastram  and 
Tessmer,  1994;  Wang  and  Schuster,  1996;  De  Lilia,  1997)  combined  with  Robertsson’s  free  surface 


48 


algorithm.  This  would  allow  the  fine  gridding  to  be  restricted  to  the  near-surface  region,  while 

coarser  gridding  could  be  used  in  other  parts  of  the  model. 

Here,  we  show  an  alternative  method  to  incorporate  surface  topography.  An  lrregular-tnangular- 
grid  method  is  presented  in  Nolte  et  al.  (1996)  which  is  based  on  a  finite-volume  approximation 
(e.g.,  Vinokur,  1989)  to  the  wave  equation.  This  method  allows  the  grid  itself  to  be  distorted  so 
that  grid  points  can  be  positioned  directly  on  interfaces.  Also,  it  deals  with  an  irregular  free  surface 
the  same  way  as  with  an  irregular  solid-solid  interface,  which  is  an  advantage  over  the  technique  of 
Muir  et  al  (1992). 

We  demonstrate  the  usefulness  of  this  method  with  numerical  examples.  First,  we  compare  the 
irregular-grid  technique  for  a  model  containing  a  sloping  interface  with  a  regular  finite-difference 
code,  and  show  that  the  former  can  model  this  case  more  accurately  than  the  latter.  We  then 
present  a  numerical  example  that  shows  how  the  technique  can  be  applied  to  a  model  with  an 
irregular  free  surface.  Finally,  using  a  model  of  realistic  surface  topography,  we  show  how  the 
topography  can  effect  regional  seismograms. 

The  Modeling  Technique 

We  have  described  the  modeling  technique  elsewhere 

(Nolte  et  al,  1996)  and  will  thus  only  review  the  essentials.  As  there  was  an  error  in  the  P- 
SV  formulation  in  the  previous  paper,  we  give  the  correct  P-SV  equations  m  the  appendix.  The 
method  essentially  amounts  to  a  discrete  approximation  to  the  wave  equation  m  integral  form,  while 
the  constitutive  relation  is  approximated  in  finite-difference  form.  The  model  is  parameterized  on 
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a  triangular  Delaunay  grid  and  its  dual  Dirichlet  grid  (see  the  appendix  and  Figure  A-2  for  the 
definition  of  these  terms).  As  the  Delaunay  grid  can  be  irregular,  grid  points  can  be  placed  directly 
on  an  interface.  In  our  irregular-grid  technique  the  medium  parameters  (A,  p  and  p)  are  defined  to 
be  constant  within  each  triangular  grid  cell  (see  Figure  A-2).  Thus,  if  grid  points  are  positioned 
exactly  on  an  interface,  the  two  media  will  be  separated  exactly  by  the  interface.  In  contrast, 
finite-difference  implementations  usually  approximate  the  model  parameters  in  an  area  around  the 
grid  point  itself,  which  in  staggered-grid  techniques  does  not  clearly  define  the  exact  position  of  an 
interface. 

As  stated  above,  the  free  surface  is  treated  the  same  way  as  any  other  interface.  On  all  other 
model  boundaries,  we  use  Higdon’s  (1986,  1987)  absorbing  boundary  condition. 


Comparison  With  Finite  Difference  for  a  Dipping  Interface 

In  order  to  test  the  performance  of  the  irregular-grid  code,  we  have  designed  the  following  exper¬ 
iment.  First,  we  model  a  case  of  a  flat  interface  (Figure  la).  For  this  type  of  model  no  stair-step 
approximation  needs  to  be  made,  therefore  a  regular  finite-difference  technique  should  give  ac¬ 
curate  results.  The  depth  to  the  interface  is  800  m,  and  the  source  is  located  120  m  above  the 
interface.  The  P  and  S  velocities  and  densities  of  the  two  layers  are  given  in  the  figure.  A  snapshot 
computed  with  finite  differences  for  the  model  in  Figure  la  is  shown  in  Figure  2a.  The  snapshot 
shows  the  horizontal  component  of  particle  velocity  (i.e.,  the  component  parallel  to  the  interface) 
for  a  point- force  source  oriented  vertically  (perpendicular  to  the  interface).  The  source  signal  is  a 
Ricker  wavelet  with  a  center  frequency  of  20  Hz.  The  snapshot  is  taken  at  a  time  of  0.2  s.  Next, 
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we  rotate  the  whole  experiment  by  30°  (Figure  lb),  i.e.,  we  rotate  the  interface,  and  the  source 
orientation.  The  source  is  therefore  still  oriented  perpendicular  to  the  interface.  For  this  model  we 
now  compute  snapshots  of  the  particle  velocity  component  parallel  to  the  interface.  Thus,  with  the 
exception  of  the  effects  from  the  free  surface  and  (imperfectly)  absorbing  boundaries,  the  snapshots 
computed  for  the  tilted  model  should  look  like  those  computed  for  the  flat  model,  but  rotated  by 
30°.  Figure  2b  shows  the  result  computed  with  the  finite-difference  scheme,  while  Figure  3b  shows 
the  irregular-grid  result.  In  the  irregular-grid  computation,  the  grid  has  been  distorted  in  such  a 
way  that  particle-velocity  grid  points  lie  exactly  on  the  dipping  interfaces.  Figure  3a  is  the  same 
as  Figure  2a.  It  is  shown  here  again  in  order  to  facilitate  the  comparison. 

Comparison  of  Figure  2b  with  Figure  2a  shows  that  while  the  overall  appearance  of  the  snapshots 
is  still  similar,  some  noticeable  differences  are  observed  particularly  in  the  region  (marked  by 
a  circle)  in  Figure  2b  where  several  waves  (the  direct  S  wave,  the  reflected  S  wave,  and  the  P 
to  S  converted  wave)  interfere,  and  where  a  critical  angle  for  the  reflected  S  is  present.  The 
amplitude- versus-angle  behavior  of  the  waves  reflected  at  the  sloping  interface  is  no  longer  modeled 
accurately.  Comparison  of  Figure  3a  with  Figure  3b,  on  the  other  hand,  shows  that  the  accuracy  of 
the  amplitude- versus-angle  behavior  for  the  tilted  interface  has  improved  considerably.  Figure  3b 
is  more  or  less  a  tilted  version  of  Figure  3a. 

This  comparison  demonstrates  that,  as  expected,  the  irregular-grid  method  can  model  irregular 
interfaces  more  accurately  than  a  standard  finite-difference  scheme. 
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A  Simple  Free-Surface  Example 


As  the  irregular-grid  technique  is  able  to  handle  irregular  solid-solid  interfaces,  it  can  be  applied 
with  confidence  to  an  irregular  free  surface.  This  is  because,  as  mentioned  above,  the  method  treats 
the  free  surface  the  same  way  as  any  other  interface.  To  demonstrate  this  capability,  we  now  apply 
it  to  a  model  with  irregular  surface  topography,  where  the  grid  is  computed  in  such  a  way  that  it 
lines  up  with  the  surface. 

First,  Figure  4  shows  an  example  of  such  a  grid  approximating  a  sloping  surface.  Note  that  only 
triangles  with  at  least  one  of  their  vertices  at  the  free  surface  are  distorted.  All  other  triangles  are 
equilateral.  For  equilateral  triangles  our  modeling  scheme  is  of  second-order  accuracy.  Otherwise 
the  accuracy  is  first-order.  For  the  gridding  used  in  Figure  4  and  in  the  example  below  the  accuracy 
is  therefore  second-order  everywhere,  except  on  the  free  surface  where  it  is  first-order.  Also  note 
that  the  grid  has  been  computed  in  such  a  way  that  none  of  the  distorted  triangles  have  angles 
greater  than  90°.  This  is  desirable,  as  it  facilitates  computation  of  the  area  elements  A  and  the 
polygons  dA  (see  the  appendix  for  the  definition  of  these  terms).  The  reason  for  this  is  that  for 
triangles  with  angles  <  90°  the  circumcenter  (the  center  of  the  circle  through  the  vertices,  see 
Figure  A-2)  is  located  within  that  triangle.  A  triangle  with  an  obtuse  angle,  on  the  other  hand, 
would  have  its  circumcircle  outside  the  triangle,  which  can  be  easily  visualized  from  Figure  A-2. 

A  model  with  irregular  surface  topography  is  shown  in  Figure  5.  The  horizontal  length  of  the 
model  is  2  km,  the  average  vertical  length  is  1.1  km.  The  topographic  relief  is  200  m.  A  single 
reflecting  interface  is  located  at  an  average  depth  of  800  m  (300  m  from  the  bottom).  The  P  and 
S  velocities  and  the  density  are  given  in  the  figure. 


Figure  6  shows  a  series  of  snapshot  taken  at  time  intervals  of  0.05  s,  starting  at  a  time  of  0.1  s. 
In  this  case,  the  source  was  a  vertically  oriented  point  force,  and  the  snapshots  show  the  vertical 
component  of  particle  velocity.  The  most  prominent  effect  of  the  irregular  surface  topography  is  the 
presence  of  scattered  surface  waves,  primarily  from  the  peak  at  about  one-third  of  the  horizontal 
model  length  from  the  left.  These  scattered  waves  can  be  seen  best  at  the  later  times.  As  can  be 
seen,  both  incident  P  waves  and  S  waves  cause  these  scattered  waves. 

Modeling  Rough-Topography  Effects  on  Regional  Wave  Propaga¬ 
tion 

We  now  use  our  modeling  technique  to  study  the  effect  of  irregular  surface  topography  on  re¬ 
gional  seismograms.  The  topography  in  our  model  (Figure  7a)  is  a  realistic  example  of  the  to¬ 
pography  that  can  be  encountered  along  profiles  in  the  Middle  East.  It  was  obtained  with  the 
online  profile  maker  from  Cornell  University’s  Middle  East  and  North  Africa  Project  database 
(http://atlas.geo.cornell.edu)  for  a  profile  in  northern  Iran.  Figure  7b  shows  the  entire  model  used 
for  our  investigation.  It  consists  of  three  layers  with  P  and  S  velocities  given  in  the  figure. 

We  computed  synthetic  seismograms  for  this  model  with  and  without  including  the  surface 
topography.  We  use  a  source  signal  with  a  center  frequency  of  0.6  Hz  for  this  computation.  The 
results  are  shown  in  Figure  8.  There  are  some  distinct  differences  in  the  waveforms  obtained  for  the 
two  cases.  In  this  example,  the  most  prominent  differences  are  on  the  vertical  component  at  larger 
arrival  times  (greater  than  70  s).  More  arrivals  with  higher  amplitudes  are  present  in  the  case  of 
surface  topography.  These  arrivals  are  most  likely  due  to  surface  waves  that  arise  from  body-wave 
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scattering  at  the  rough  surface,  the  same  effect  we  observed  in  the  irregular-surface  example  in  the 
previous  section. 

Our  comparison  demonstrates  that  a  pronounced  surface  topography,  such  as  is  present  in  many 
regions,  e.g.  the  Middle  East,  will  have  a  distinct  effect  on  observed  waveforms,  and  should  thus 
be  included  in  the  modeling.  Our  newly  developed  method  provides  a  valuable  tool  for  this  task. 


Discussion 

The  memory  requirements  of  the  technique  have  not  been  addressed  thus  far.  The  irregular  grid 
requires  more  memory  per  grid  point  then  a  regular-grid  method.  The  additional  quantities  that 
need  to  be  stored  for  each  grid  point  are  its  coordinates  and  the  indices  of  its  neighbors.  Also,  the 
total  number  of  stress  points  will  be  larger  than  the  total  number  of  velocity  points  (the  exact  ratio 
of  these  two  quantities  depends  on  the  grid  irregularity).  On  the  other  hand,  the  method  allows 
variable  grid-point  spacing,  so  that  the  spacing  can  be  locally  scaled  to  the  medium  velocity  as  well 
as  to  the  amount  of  detail  with  which  one  desires  to  model  particular  areas  of  a  model.  For  example, 
for  a  model  in  which  regions  with  strongly  different  velocities  are  present,  the  irregular-grid  method 
would  require  fewer  grid  points  than  a  second-order  finite-difference  program  and  would  thus  allow 
the  model  size  to  be  larger. 

A  potential  extension  of  the  code  to  3-D  is  easy  in  principle.  Taking  the  method  to  3-D  would 
require  tetrahedral  instead  of  triangular  Delaunay  cells.  In  practice,  one  of  the  main  difficulties  is 
the  problem  of  the  grid  computation  itself,  i.e.,  fitting  a  tetrahedral  grid  to  a  model  in  such  a  way 
that  it  lines  up  with  the  interfaces  and  that  the  grid  spacing  is  scaled  to  the  medium  velocities. 
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However,  the  development  of  efficient  3-D  gridding  algorithms  is  an  area  of  extensive  research  (e.g., 
the  GOCAD  Research  Program),  so  that  existing  software  can  be  used  to  perform  this  task. 

Conclusions 

In  this  study,  we  have  demonstrated  the  ability  of  the  irregular-grid  technique  to  accurately  model 
irregular  interfaces.  A  comparison  with  a  finite-difference  technique  showed  that  the  irregular-grid 
method  produced  superior  results,  as  no  stair-step  approximation  to  interfaces  is  made.  We  also 
showed  the  applicability  of  the  method  to  an  irregular  free  surface.  A  numerical  test  showed  that 
surface  irregularities  gave  rise  to  numerous  scattered  waves.  We  observed  the  same  scattering  effect 
in  regional  synthetics  computed  for  a  model  with  realistic  surface  topography. 
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Appendix 


A  more  detailed  description  of  the  modeling  technique  is  given  in  Nolte  et  al.  (1996).  Here,  we  only 
give  a  brief  review  of  the  method  for  the  P-SV  case  and  correct  some  errors  of  the  earlier  paper. 

Figure  A-l  defines  some  terms  used  in  the  equations  below.  For  an  area  A  with  boundary  d A, 
we  introduce  local  coordinates  n  and  s,  normal  and  tangential  to  dA,  respectively.  With  <f>  denoting 
the  angle  between  n  and  the  Cartesian  x  direction,  the  local  and  Cartesian  coordinates  are  related 
through: 

n  =  x  cos  (j)  +  z  sin  (j> 
s  =  —  arsing  +  zcos<t> 
x  =  n  cos  <j>  —  s  sin  4> 

z  =  nsm<f>  +  scos<f>  (1) 


Using  the  above  definitions  we  write  the  wave  equation  in  integral  form  as 


( b  ds  (ann  cos  <j>  -  i 7ns  sin  <f>) 
J  dA 

<p  ds  (ann  sin  <f>  +  ans  cos  (j)) . 
J  dA 


(2) 


Here,  vx  and  vz  are  the  x  and  z  components,  respectively,  of  the  particle  velocity,  while  ann  and 
<rns  are  the  normal  and  tangential  stresses  on  dA. 

The  constitutive  relation  (Hooke’s  law)  can  be  written  as 


dcnn 

dt 

d&ns 

dt 


ds 


,dvx  dvz  .  ,,dvz  dvx  .  . 

(A  +  2m)(-^  C0S(f>+~d^  sm  ^  +  x(~d7 cos  ^  “  sin  w 

,dvz  dvx.  ,  ,dvz  dvx  . 

(_,  +  _Z)c°s0  +  (_-_)sm^ 


(3) 
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In  order  to  discretize  these  equations,  the  model  is  parameterized  using  Delaunay  and  Dirichlet 
grids  (Delaunay,  1934;  Dirichlet,  1850)  as  illustrated  in  Figure  A-2.  A  Delaunay  grid  is  a  triangular 
grid  that  connects  nearest  neighbors.  This  nearest-neighbor  condition  means  that  a  circle  through 
the  three  vertices  of  any  triangle  (a  so-called  circumcircle)  does  not  enclose  any  other  vertex. 
The  corresponding  Dirichlet  grid  is  obtained  by  connecting  the  centers  of  the  circumcircles  (the 
circumcenters )  of  neighboring  Delaunay  triangles. 

We  now  specify  particle  velocity  at  the  Delaunay  vertices,  and  use  the  area  enclosed  by  the 
polygon  around  any  particle-velocity  point  as  the  area  A  of  Figure  A-l  while  we  use  the  polygon 
itself  as  its  boundary  dA.  We  assume  that  the  medium  parameters  (A,  p,  and  p)  as  well  as  all 
stresses  are  constant  within  each  Delaunay  triangle. 

The  notations  used  are  shown  in  Figure  A-3a  and  A-3b.  In  Figure  A-3b,  we  have  introduced 
the  velocities  vi,V2,v$,  and  V3,  for  convenience.  As  can  be  seen,  V3  is  the  velocity  vy.fe-i],  and  v2  is 
V[7,*+i];  vl  is  interpolated  from  Vj,  v\jtk}i  and  v2\  V4  is  interpolated  from  Vj,  V3,  and  Denoting 

the  time  index  by  m  we  obtain  the  discrete  form  of  the  wave  equation  as 


m  _ 


vx\j] 


vm 


m— 1  , 

V*\j)  + 


\  Nk 

■f-r )  ■ -  «21  ^ 

1  A\jMP\iM  J  fc=i 


Efi 


+l\j4a™VM  cos  <j>-  aJyM  sin<j>}) 

\  Nk 


TO-i  , 
+ 


At 


+lv4a™m/2]  Sin  ^  +  anslM/2]  C0S  ^)- 


sin  *  +  an}m/2]  cos  $ 


(4) 
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We  obtain  the  discrete  form  of  the  constitutive  relation  as 


.+  [m+1/2] 

_  +  [m— 1/2] 

nn  \j,k] 

nn  [?,fc] 

-  [m-f  1/2] 

_  -  [m-l/2] 

nn  [j,k] 

y  nn  [7,/c] 

-f  [m+1/2] 

__  _+  [m_1/2] 

ns  [j,k] 

°ns  [j,A:] 

[m+1/2] 

_  —  [m— i/2] 

ns  [j,fc] 

ns  [?,/c] 

+ 


At  (A+  +  2M+)(Aun)b-fc]  +  A+4^(A vs)+ 


dbM 


em 


+  -y—  (A  +  2/i  )(Avn)b>fc]  +  A  +  -3— (Ava)biJb] 

ey,fc]  '  ' 

+ ■£r,l,+iAv')M + ^-fi+(A,,")S.‘i 

aM  %M 

At  _/A  .  At 

+  —^  (A^)[i)fc]  +  — M  (Avn)m 


lm 


(5) 


The  indices  +  and  -  refer  to  the  two  different  grid  cells,  respectively,  bordering  on  the  line 
from  point  j  to  neighbor  k  (see  Figure  A-3b).  In  Eq.  (5)  we  have  used  the  following  abbreviations 
for  the  normal  velocity  derivatives 

Avn  =  (Ck  -  v™j)  cos  4>  4-  (Cfc  —  v™j)  sin  <j>] 

Avs  =  -(Ck  -  v™)  sm4>  +  (Ck  -  C;)  cos  <j>\.  (6) 

and  for  tangential  velocity  derivatives  (see  Figure  A-2b) 

(Avn)+  =  (C2  -  Ci)  cos  4>  +  (v™2  -  Cl)  sin  <f>] 

(Avn)~  =  (Cl  ~  Cs) cos  ^  +  (Ct  -  Cs) sin 
(Avs)+  =  -(C2  “  Ci) sin <t>  +  (C  ~  Ci) cos # 

(Avs)“  =  -(C4-C3)sm^  +  (C4-C3)cos#  (7) 

As  can  be  seen  in  the  above  equations  the  particle  velocities  are  updated  at  time  steps  m  =  1, 
2,  ...  whereas  stresses  are  updated  at  time  steps  m  =  1/2,  3/2,  . . .  analogously  to  a  staggered-grid 
finite-difference  approach. 
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Figure  1:  Models  used  for  comparing  the  irregular-grid  method  with  finite  differences.  The  P  and 
S  velocities  (in  kms-1)  and  densities  (in  gem-3}  are  given  in  the  figure,  (a)  A  fiat  interface 
is  located  at  a  depth  of  800  m.  The  source  is  a  point  force  located  120  m  above  the  interface. 
The  source  orientation  is  perpendicular  to  the  interface  as  indicated  by  an  arrow'  in  the  figure, 
(b)  The  experiment  is  rotated  by  30°.  The  source  orientation  is  again  indicated  by  an  arrow-  in 
the  picture. 
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Figure  2:  Snapshots  computed  with  finite  differences  for  (a)  the  model  of  Figure  la  and  (b)  the 
model  of  Figure  lb.  The  snapshots  show  the  component  of  particle  velocity  parallel  to  the 
interface  at  the  time  of  0.2  s.  The  waves  identified  in  (a)  are  the  direct  P  wave  (1),  the  direct  S 
wave  (2),  the  P-to-S  converted  wave  (3),  and  the  reflected  S  wave  (4).  The  region  enclosed  by 
the  circle  in  (b)  is  inaccurately  modeled  with  finite  differences. 


Figure  3:  Same  as  Figure  2,  but  showing  the  snapshot  computed  with  the  irregular-grid  method 
(b;  for  the  tilted-mterface  model.  Note  the  improved  accuracy,  particularly  in  the  region  where 
the  finite-differences  results  were  inaccurate  (the  region  circled  in  Figure  2b). 


',A\  A  A  A  A  A  A  A  A/V  V  V  \/  \/  \/  \/  w  w  w  w  ^ 


Figure  4:  Example  of  a  triangular  grid  that  is  lined  up  with  an  irregular  free  surface.  All  triangles 
except  those  that  touch  the  free  surface  are  equilateral.  No  angles  greater  than  90°  axe  present 
anywhere  in  the  grid. 
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Figure  5:  Model  with  irregular  topography.  The  P  and  S  velocities  (in  kxns-1)  and  densities  (in 
g  cm-3)  are  given  in  the  figure.  The  length  of  the  model  is  2  km,  its  height  varies  between  1 
km  and  1.2  km.  This  model  is  the  input  for  the  snapshots  of  Figure  6. 
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Figure  6:  Snapshots  computed  for  the  model  in  Figure  5.  The  snapshots  show  the  vertical  compo¬ 


nent  of  particle  velocity  for  at  time  intervals  of  0.05  s,  starting  at  0.1  s. 
point  force. 


The  source  is  a  vertical 
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Figure  7:  Irregular  surface  topography  (a)  and  three-layer  model  (b)  for  regional  wave  propagation 


simulation.  The  source-receiver  offset  is  259  km.  The  source  depth  is  10  km. 


No  Topography,  Horizontal  Component 


No  Topography,  Vertical  Component 


Topography,  Horizontal  Component 


Time  (s) 


Figure  8:  Synthetic  waveforms  computed  for  the  model  in  Figure  7  with  and  without  accounting 
for  surface  topography. 


71 


Figure  A-l:  Definition  of  the  area  A,  its  boundary  dA,  the  local  coordinates  n  and  s,  and  the  angle 
<fr  used  in  the  equations  in  the  appendix. 
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Figure  A-2:  Schematic  of  a  Delaunay  grid  and  its  dual  Dirichlet  grid.  The  circumcircle  through  the 
vertices  of  any  Delaunay  triangle  does  not  contain  any  other  triangle  vertices.  The  circumcenters 
are  the  vertices  of  the  Dirichlet  grid. 
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Figure  A-3:  Sketch  of  the  terms  used  in  the  discrete  equations,  (a)  A  particle-velocity  point  j 
with  six  neighbors  k  is  shown.  Stress  is  integrated  over  the  dashed  polygon.  Ayik]  is  an  area 
element,  and  dy  kj  is  the  distance  from  point  j  to  neighbor  k.  (b)  Line  segments  ly>k]  over  which 
stress  is  summed;  velocity  points  vi,V2,t>3,  and  v.\  and  lengths  e\j,k]  used  for  the  computation  of 
tangential  velocity  differences. 
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LAMONT-DOHERTY  EARTH  OBSERVATORY 
PALISADES,  NY  10964 


DAVID  RUSSELL 
HQ  AFTAC/TTR 
1030  SOUTH  HIGHWAY  A1 A 
PATRICK  AFB,  FL  32925-3002 


SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5704 

MS  0979,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0979 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5704 

MS  0655,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0655 

THOMAS  SERENOJR. 

SCIENCE  APPLICATIONS  INTERNATIONAL 

CORPORATION 

10260  CAMPUS  POINT  DRIVE 

SAN  DIEGO,  CA  92121 

ROBERT  SHUMWAY 
410  MRAK  HALL 
DIVISION  OF  STATISTICS 
UNIVERSITY  OF  CALIFORNIA 
DAVIS,  CA  956 16-8671 


JEFFRY  STEVENS 
MAXWELL  TECHNOLOGIES 
8888  BALBOA  AVE. 

SAN  DIEGO,  CA  92123-1506 


TACTEC 

BATTELLE  MEMORIAL  INSTITUTE 
505  KING  AVENUE 

COLUMBUS,  OH  43201  (FINAL  REPORT) 


DELAINE  REITER 
SENCOM  CORP. 

73  STANDISH  ROAD 
WATERTOWN,  MA  02172 


MICHAEL  RITZWOLLER 
DEPARTMENT  OF  PHYSICS 
UNIVERSITY  OF  COLORADO 
CAMPUS  BOX  390 
BOULDER,  CO  80309-0309 


CHANDAN  SAIKIA 

WOOODWARD-CLYDE  FEDERAL  SERVICES 
566  EL  DORADO  ST.,  SUITE  100 
PASADENA,  CA  91 101-2560 


SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  93 11 

MS  1 159,  PO  BOX  5800 
ALBUQUERQUE,  NM  87 1 85- 1 1 59 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5736 

MS  0655,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0655 


AVI  SHAPIRA 
SEISMOLOGY  DIVISION 

THE  INSTITUTE  FOR  PETROLEUM  RESEARCH  AND 
GEOPHYSICS 

P.O.B.  2286  NOLON  58122  ISRAEL 


DAVID  SIMPSON 
IRIS 

1200  NEW  YORK  AVE.,  NW 
SUITE  800 

WASHINGTON  DC  20005 


BRIAN  SULLIVAN 
BOSTON  COLLEGE 
INSITUTE  FOR  SPACE  RESEARCH 
140  COMMONWEALTH  AVENUE 
CHESTNUT  HILL,  MA  02167 


NAFI  TOKSOZ 

EARTH  RESOURCES  LABORATORY,  M.I.T. 
42  CARLTON  STREET,  E34-440 
CAMBRIDGE,  MA  02142 


LAWRENCE  TURNBULL 

ACIS 

DCI/ACIS 

WASHINGTON  DC  20505 


GREG  VAN  DER  VINK 
IRIS 

1200  NEW  YORK  AVE.,  NW 
SUITE  800 

WASHINGTON  DC  20005 
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FRANK  VERNON 

UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 
SCRIPPS  INSTITUTION  OF  OCEANOGRAPHY 
9500  GILMAN  DRIVE 
LA  JOLLA,  CA  92093-0225 


TERRY  WALLACE 
UNIVERSITY  OF  ARIZONA 
*,  0225  DEPARTMENT  OF  GEOSCIENCES 

BUILDING  #77 
TUCSON,  AZ  85721 


JILL  WARREN 

LOS  ALAMOS  NATIONAL  LABORATORY 
GROUP  NIS-8 
P.O.  BOX  1663 

LOS  ALAMOS,  NM  87545  (5  COPIES) 


DANIEL  WEILL 
NSF 

EAR-785 

4201  WILSON  BLVD.,  ROOM  785 
ARLINGTON,  VA  22230 


RUSHAN  WU 

UNIVERSITY  OF  CALIFORNIA  SANTA  CRUZ 
EARTH  SCIENCES  DEPT. 

1156  HIGH  STREET 
SANTA  CRUZ,  CA  95064 


JIAKANG  XIE 
COLUMBIA  UNIVERSITY 
LAMONT  DOHERTY  EARTH  OBSERVATORY 
ROUTE  9W 

PALISADES,  NY  10964 


JAMES  E.  ZOLLWEG 
BOISE  STATE  UNIVERSITY 
GEOSCIENCES  DEPT. 

1910  UNIVERSITY  DRIVE 
BOISE,  ID  83725 
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